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ABSTRACT 

We use a new multiannulus planetesimal accretion code to investigate the evolu- 
tion of a planetesimal disk following a moderately close encounter with a passing star. 
The calculations include fragmentation, gas and Poynting-Robertson drag, and veloc- 
ity evolution from dynamical friction and viscous stirring. We assume that the stellar 
encounter increases planetesimal velocities to the shattering velocity, initiating a colli- 
sional cascade in the disk. During the early stages of our calculations, erosive collisions 
damp particle velocities and produce substantial amounts of dust. For a wide range 
of initial conditions and input parameters, the time evolution of the dust luminosity 
follows a simple relation, L^/L^ = Lo/[a + (t/td)^]- The maximum dust luminosity Lq 
and the damping time t^ depend on the disk mass, with Lq oc and t^ oc Mj^. For 
disks with dust masses of 1% to 100% of the 'minimum mass solar nebula' (1-100 M<^ at 
30-150 AU), our calculations yield td ~ 1-10 Myr, a ^ 1-2, P = 1, and dust luminosi- 
ties similar to the range observed in known 'debris disk' systems, Lq ~ 10^^ to 10^^. 
Less massive disks produce smaller dust luminosities and damp on longer timescales. 
Because encounters with field stars are rare, these results imply that moderately close 
stellar flybys cannot explain collisional cascades in debris disk systems with stellar ages 
of ~ 100 Myr or longer. 

Subject headings: planetary systems - solar system: formation - stars: formation - 
circumstellar matter 



1. INTRODUCTION 



Many nearby main sequence stars have thermal emission from cold dust (Backman &: Paresce 
1993; Artymowicz 1997; Lagrange et al. 2000). These 'debris disk' systems often have a disk-like 
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or ring morphology at optical or infrared wavelengths, with typical radii of 50-1000 AU (e.g., Smith 

& Terrile 1984; Greaves ct al. 1998; Holland et al. 1998; Jayawardhana ct al. 1998; Kocrner 
et al. 1998; Schneider et al. 1999; Weinberger et al. 1999). The near-IR and far-IR excess 
emission is a small fraction of the stellar luminosity; the ratio of the total dust luminosity to the 
stellar luminosity is L^/L^ ~ 10"^ to 10~^ in most systems (Aumann et al. 1984; Fajardo-Acosta 
et al. 1999; Habing et al. 1999; Spangler et al. 2001). Recent studies indicate that the dust 
luminosity declines on timcscalcs of 10-100 Myr (Barrado y Navascues et al. 1999; Habing et al. 
1999; Song et al. 2000b). Kalas (1998) and Spangler et al. (2001) propose power-law relations 
for this decline, L^i oc (Kalas 1998) and L^i oc (Spangler et al. 2001). 

In current theories, dynamic processes produce dust emission in a debris disk. If the disk 
contains small dust grains with sizes of 1-100 ;um and a total mass in small grains of ~ 0.01 Mq^, 
radiative transfer models can explain both the scattered-light images and the spectral energy dis- 
tributions of most systems (e.g., Backman k. Paresce 1993; Artymowicz 1997, Greaves, Mannings, 
& Holland 2000b). However, radiation pressure and Poynting-Robertson drag remove small grains 
from the disk on short timescales, ^ 1-10 Myr, compared to the age of a typical system. Collisions 
between larger bodies can replenish the small grain population when their relative velocities are 
large enough to begin a 'collisional cascade,' where 1-10 km planetesimals are ground down into 
smaller and smaller bodies. If planetesimal velocities are close to the 'shattering velocity' of ~ 200 
m s"-*^ and if the mass in 1-10 km planetesimals is ~ 100 M®, collisional cascades can maintain 
the observed small grain population over timescales of ~ 100 Myr or longer (Habing et al. 1999; 
Greaves, et al. 2000b; Kenyon & Bromley 2001). 

At least two mechanisms can produce large collision velocities in a debris disk. During the 
early evolution of the disk, there is a strong coupling between solid bodies and the gas; circular- 
ization is efficient and collision velocities are low. As solid bodies in the disk merge and grow, 
impulsive encounters with passing stars or stirring by large planets embedded in the disk can in- 
crease collision velocities (Artymowicz et al. 1989; Mouillet et al. 1997; Ida et al. 2000; Kalas 
et al. 2000). Kenyon k. Bromley (2001) show that embedded planets with radii of 500 km or 
larger can dynamically heat smaller bodies to the shattering velocity on timescales of 10-100 Myr. 
Collisions between 1-10 km objects can then replenish the small grain population for ~ 500 Myr 
or longer (sec also Kenyon 2002). Ida et al. (2000) and Kalas et al. (2000) use n-body simulations 
to demonstrate that close encounters between the disk and a passing star can excite large particle 
velocities in the disk. Because the growth of large planets in the outer disk is unlikely, these 'stellar 
fiybys' are an attractive way to induce collisional cascades in the large disks observed in (3 Pic and 
other nearby stars. The evolution of a particle disk after the stellar flyby is uncertain. Collisions 
of small bodies in the disk may damp planetesimal velocities and thus halt the collisional cascade. 
Kenyon &: Luu (1999) derive damping times of 1-10 Myr for 1-100 m planetesimals in the Kuiper 
Belt. If this timescale is typical, then stellar fiybys can explain debris disks only if the encounters 
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axe frequent. 

Here we investigate the evolution of a planetesimal disk following a moderately close stellar 
flyby. We use a new multi-annulus planetesimal evolution code to compute the coUisional and 
dynamical evolution of small bodies in an extended disk surrounding an intermediate mass star. 
Calculations for disks with dust masses comparable to the 'minimum mass solar nebula' (~ 100 
of solid material at 30-150 AU) show that two body collisions can damp planetesimal velocities on 
timescales of 1-10 Myr. Less massive disks damp on longer timescales. Small dust grains generated 
in the collisional cascade produce a significant luminosity, with L^/L^ ~ 5 x 10~^ for a minimum 
mass solar nebula. Less massive disks yield smaller dust luminosities. These results suggest that 
stellar flybys cannot initiate luminous, long-lived collisional cascades in a planetesimal disk. Unless 
all known debris disk systems have had a stellar encounter in the past 10-100 Myr, a nearby binary 
companion star or planets embedded in the disk probably initiated the collisional cascades in these 
systems. 

We outline the model in §2, describe the calculations in §3, and conclude with a brief discussion 
in S4. 



2. THE MODEL 

We treat planetesimals as a statistical ensemble of masses with a distribution of horizontal and 
vertical velocities about a Keplerian orbit (Safronov 1969; Wetherill & Stewart 1989; Kenyon & 
Luu 1998, 1999, and references therein). The model grid contains N concentric annuli centered 

at heliocentric distances o,;. Each annulus has an inner radius at Oj — 5ai/2 and an outer radius at 
tti + 5ai/2. The midpoint of the model grid is at a heliocentric distance Umid- Calculations begin 
with a differential mass distribution n{mik) of bodies with horizontal and vertical velocities hik{t) 
and Vik{t). 

To evolve the mass and velocity distributions in time, we solve the coagulation and the Fokker- 
Planck equations for an ensemble of masses undergoing inelastic collisions, drag forces, and dynam- 
ical friction and viscous stirring. We approximate collision cross-sections using the particle-in-a-box 
formalism and treat gas and Poynting-Robertson drag with simple analytic formulae (see the Ap- 
pendix). Kenyon & Bromley (2001) outline our treatment of elastic gravitational interactions (see 
also Spaute et al. 1991; Weidenschilling et al. 1997). 

We allow four collision outcomes (Greenberg et al. 1978; Davis et al. 1985; Barge k, Pellat 
1993; Wetherih & Stewart 1993): 

1. Mergers - two bodies collide and merge into a single object with no debris; 

2. Cratering - two bodies merge into a single object but produce debris with mass of 20% or less 
of the mass of the merged object; 
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3. Rebounds - two bodies collide and produce some debris but do not merge into a single object; 
and 

4. Disruption - two bodies collide and produce debris with a mass comparable to the mass of the 
two initial bodies. 

The appendix of Kenyon & Luu (1999) describes the algorithm for specifying the outcome as a 
function of the collision velocity, the tensile strength of the colliding bodies Sq, and other physical 
parameters. For the calculations described in this paper, all collisions result in debris; rebounds 
are unimportant. For most calculations, we adopt Sq = 2 x 10^ erg g"^ and a crushing energy 
Qc = 5 X l{f erg g~^, appropriate for icy objects at large distances from the central star (see 
Kenyon &; Luu 1999, and references therein). 

The initial conditions for these calculations are appropriate for a disk with an age of ~ 10 Myr. 

Wc consider systems of N annuli in disks with dmid — 

35-140 AU and 5ai/ai = 0.01-0.025. The 
disk is composed of small planetesimals with radii of ~ 1-100 m (see below). The particles have 
an initial mass distribution ni{mk) in each annulus and begin with eccentricity eo and inclination 
io = eo/2. Most of our models have cq independent of af, some models have a constant initial 
horizontal velocity in each annulus. We consider moderately strong perturbations with cq = 0.01- 
0.06. Very close stellar flybys can produce cq > 0.1 (e.g., Ida ct al. 2000; Kalas ct al. 2000; 
Kobayashi & Ida 2001); our adopted cq is 10-100 times larger than for a cold planctesimal disk, 
where eo < lO"-'^ (Kenyon & Lmi 1999). Wc assume a power law variation of the initial surface 
density of solid material with heliocentric distance, Sj = So(ai/ao)^"; models with n = 1.5 are 
'standard'. 

We assume a stellar flyby instantaneously raises the orbital eccentricity and inclination of 
all planetesimals in the disk. Although not strictly accurate, this approximation is valid if the 
encounter between the passing star and the disk is short compared to the damping time. Any 

passing star that is not bound to the central star of the disk satisfies this requirement (Larwood 
1997; Mouillet et al. 1997). A bound companion star with an orbital semi-major axis of 1000 AU 
or less can excite large planctesimal velocities over long timescales and may counteract damping. 
Embedded planets with radii of 500 km or larger can also counteract damping (Kenyon & Bromley 
2001). We plan to consider both situations in future papers. 

Our calculations follow the time evolution of the mass and velocity distributions of objects 
with a range of radii, rj^ = rmin to r^fc = Vmax- The upper limit rmax is always larger than the 
largest object in each annulus. To save computer time, we do not consider small objects which 
do not affect significantly the dynamics and growth of larger objects. Usually this limit is rmin = 
10 cm. Erosive collisions produce objects with r^fc < r„iin- These particles are 'lost' to the model 
grid. For most conditions, lost objects are more likely to be ground down into smaller objects than 
to collide with larger objects in the grid. Thus, our neglect of lost particles with rik < Vmin is a 
reasonable approximation. 
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We do not consider planetesimals in the inner disk. Test calculations show that collisions 
damp particle velocities on short timescales, < 10^ yr at aj < 5 AU. For a stellar flyby with an 
encounter velocity of 1-10 km s~^, these timescales are comparable to the duration of the flyby 
(see Larwood 1997; Mouillet et al. 1997). Our approximation of an instantaneous impulse to the 
disk in then invalid. The short damping timescales for the inner disk also suggest that the chances 
of finding a disk in this state are small. We therefore consider the evolution of the outer disk in 
these calculations. 

For a cold disk of planetesimals with cq < 10~^, the initial conditions of our calculations 
are sufficient to produce 1000 km or larger planets at 30 AU on timescales of ~ 10 Myr for a 
minimum mass solar nebula (e.g., Stern &; Colwell 1997a; Kenyon & Luu 1999; Kenyon et al. 
1999). When a stellar flyby raises particle eccentricities by factors of 10-100, collisions between 
planetesimals promote erosion instead of growth by mergers. Our goal is to identify perturbations 
where collisional damping can reduce particle eccentricities and promote growth before erosive 
collisions reduce planetesimals to dust. 



3. CALCULATIONS 

3.1. Cascades in Rings 

To understand the response of a planetesimal disk to a stellar flyby, we begin by considering 
a small portion of a disk surrounding a 3 Mq star. The ring consists of bodies with rmin = 10 cm 
to rmax in 16 annuli with Sai/ui = 0.01. Each body has eo = 0.02 and a mass density of 1.5 g 
cm~^. The mass spacing factor is Si = mk+i/mk = 2 for all batches; the initial mass distribution 
has equal mass per batch. The surface density of solid material is Sd(a) = 0.15 x(a/35 AU)^'^/^ 
g cm~^, where a; is a dimensionless constant. Models with x = 1-2 have a total mass in solid 
material similar to the 'minimum mass solar nebula' at 2-20 AU (Weidenschilling 1977a; Hayashi 
1981). 

We calculate the collisional evolution for a set of models with x = 10^^ to 1 and r^ax = 
10 m to 1000 km at amid = 35 AU, 70 AU, and 140 AU. For computational speed, we adopt the 
Wetherill &; Stewart (1993) fragmentation algorithm. Calculations with the Davis et al. (1985, 
1994) algorithm evolve on a faster timescale when the collisional debris receives more kinetic energy 
per unit mass than the merged object (see Kenyon & Luu 1999). Collisional debris rarely has less 
kinetic energy per unit mass than the merged object (Davis et al. 1985, 1994). The fragmentation 
parameters for our standard models are Qc = 5 x 10'' erg g^^ (crushing energy), Sq = 2 x 10^ erg 
g~^ (impact strength), and Vj = 1 m s^^ (the minimum velocity for cratering). 

Figures 1~2 outline the evolution of a model with x = 0.1 and rmax = 10 m at amid = 35 AU. 
The 16 annuli in this standard model contain an initial mass of 4 X lO^^^ g (0.67 Me). Erosive 
collisions dominate the early stages of this calculation. Collisions between small bodies with = 
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10-100 cm are completely disruptive; collisions between the largest bodies yield some growth and 
substantial cratering. Throughout the first 10^ yr of the calculation, the initial 'mass loss rate' is 
large: ~ 10^^ g yr^-*^ is lost to particles with r/. < Tmin = 10 cm. This rate falls uniformly with time 
thereafter. Disruptive collisions between the 'lost' particles probably produce micron-sized grains 
which are ejected from the system on short timescales, < 10^ yr (e.g., Krivov et al. 2000). 

As the calculation continues past 10^ yr, collisional damping reduces the velocities of the mid- 
sized particles with ~ 1 m. Collisions between these and the largest bodies in the grid produce 
growth and some debris. The debris produced from these and other collisions flattens the cumulative 
number distribution from the initial Nc oc to Nc oc m~°'^^. For the rest of the calculation, 
the power law slope of the number distribution for the small bodies remains close to —0.83 in each 
annulus. This slope is a standard result of coagulation when collisions produce debris (Dohnanyi 
1969; Wetherih & Stewart 1993; Williams k Wetherih 1994; Tanaka et al. 1996; Kenyon & 
Luu 1999). As collisional damping continues to reduce the velocities of the mid-sized bodies, the 
largest bodies grow more rapidly. By t = 1 Myr, the largest bodies have = 25 m (Figure 1, 
middle left panel). Although the largest and smallest bodies retain a large fraction of their initial 
eccentricities, the orbits of the mid-sized bodies are a factor of ~ 4 less eccentric than their initial 
orbits (Figure 1; middle right panel). As growth proceeds, dynamical friction considerably reduces 
the velocities of the largest objects. Collisions damp the velocities of the smaller bodies. By t = 
3 Myr, the largest objects have = 100 m and most objects have Sk < 0.1 cq (Figure 1, lower 
panels). From previous results, the largest bodies soon begin 'runaway growth' and reach sizes of 
10 100 km or larger on timescales of ~ 30-100 Myr (Greenberg ct al. 1978; Kenyon & Luu 1999; 
Kenyon et al. 1999). Because the production of debris - already below ~ 10^"^ g yr~^ — drops 
rapidly during runaway growth, we stopped this calculation at t = 10 Myr. 

Figure 2 shows the time evolution of two 'observables', the radial optical depth and the fraction 
of stellar luminosity intercepted and reprocessed by solid material in the grid, L/L^. We follow 
Kenyon et al. (1999) and calculate the radial optical depth of particles in the grid using the 
geometric optics limit. For particles with < 10 cm, we estimate a minimum optical depth Tmin 
assuming Poynting-Robertson drag is efficient at removing the smallest objects and a maximum 
optical depth T^ax assuming Poynting-Robertson drag removes none of the smallest objects. The 
dust luminosity follows directly from the solid angle of the dusty annulus as seen from the central 
star, L/L^ = rHd/a^ where Hd is the scale height of the dust. The appendix describes our derivation 
of the radial optical depth and the reprocessed stellar luminosity in more detail (see also Krivov et 
al. 2000). 

At t = 0, the optical depth is large because the dust production rate is large. The optical depth 
of the smallest grains, Tmin and Tmaxj fells with the declining dust production rate and drops by ~ 
two orders of magnitude by t = 10 Myr. For t < 1 Myr, the optical depth of the mid-sized bodies 
grows with time. Collisional damping reduces the scale height of these particles, which increases 
their number density and thus their radial optical depth. By t = 10 Myr, the largest bodies begin 
to grow rapidly and falls back to its initial level. The optical depth continues to drop as runaway 
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growth concentrates more and more mass in the largest bodies. 

The evolution of the dust luminosity depends on the evolution of the optical depth and the 
dust scale height^ Because collisional damping reduces particle velocities throughout the calculation, 
the dust scale height falls monotonically with time. Unless the optical depth rises substantially, 
decreasing dust scale heights result in a declining dust luminosity. The small rise in Td for t < 
1 Myr does not compensate for the decline in the scale height: the dust luminosity thus declines 
with time. At t = 0, the dust luminosity roughly equals the observed luminosities of debris disk 
systems with ring-like geometries, L/L^ ~ 10~^ (e.g., Augereau et al. 1999b; Greaves et al. 1998; 
Jayawardhana et al. 1998; Koerner et al. 1998). By the end of the calculation, our most optimistic 
estimate of the dust luminosity falls below the smallest observed luminosity in known debris disk 
systems. 

Although our predicted luminosities are close to those observed, the radial optical depths in our 
calculations, ^ 10^^ for the small particles and ~ 10~^ large particles in the grid, are significantly 
larger than values of ~ 10~^ or less typically quoted for debris disk systems. Published optical 
depths are usually derived from the ratios of the dust luminosity to the stellar luminosity and 
assume the dust is arranged in a spherical shell, L/L^, = r(l — w), where u> is the albedo. We adopt 
a more appropriate ring geometry; the dust luminosity is L/L^ = r(l — u))H/a. For rings with 
H/a ~ 0.01-0.1, the radial optical depths must be ~ 10—100 times larger through the midplane 
of a ring than through a spherical shell. With this scaling, our optical depth estimates agree with 
published estimates for debris disk systems. 

The luminosity evolution for each curve in Figure 2 is well-fit by a simple function 

L/L^ = ^^V^ , (1) 

' (x+(t/tdY ^ ' 

where Lq is the initial dust luminosity. Because collisional damping reduces the eccentricities of 
the mid-sized particles by a factor of ~ 2 at t = t^, we call td the 'damping time.' At late times, 
the luminosity follows a power law with index (3. We derive (3 = 1.35 ib 0.02 for the luminosity of 
small particles, Lmin and L^ax-, and (3 = 1.05 it 0.02 for the luminosity of large particles, L^,. 

Calculations with other initial conditions produce similar results but on different timescales. 
High velocity collisions convert 25% or more of the initial mass into small particles 'lost' to the 
grid. Cratering produces most of the mass loss; catastrophic disruption is responsible for < 10% of 
the lost mass. Collisional damping reduces the velocities of the mid-sized particles before cratering 
and catastrophic disruption convert all of the initial mass into debris. Because the collision lifetime 



^We do not attempt to model the evolution of the dust luminosity during the stellax flyby. Because we assume 
a size distribution for small particles, the luminosity at t = in our models assumes that these particles reach the 
collisional size distribution, Nc oc m^"'*^, instantaneously. In a more realsitic calculation, the timescale to achieve 
this size distribution, ~ 10'* yr to 10^ yr, is short compared to the damping time. Thus, our model provides a 
reasonably good approximation to the luminosity evolution following the flyby, but it makes no prediction for the 
change in L/Lo resulting from the flyby. 
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of a particle scales with the particle density, collisional damping is more effective in the more 
massive grids. Once the eccentricities of the mid-sized particles reach ~ 5 x 10~^, collisions begin 
to produce larger objects and less debris. Dynamical friction then reduces the velocities of the 
largest bodies, while collisional damping continues to 'cool' the mid-sized and smaller bodies. The 
most massive bodies eventually reach velocities small enough to initiate runaway growth, where the 
largest bodies grow rapidly and the mass loss rate is negligible (Kenyon & Luu 1999). We plan to 
describe runaway growth and the formation of planets in a future paper. 

Figure 3 outlines the evolution of a model with x = 0.1 and r^ax = 10 km at amid = 35 AU. 
This model has fewer small particles than a grid with Vmax = 10 m; the smaller collision rate thus 
leads to a smaller luminosity and a smaller damping rate. The damping of the smallest particles is 
complete hy t = 1 Myr, but the largest particles still retain most of their initial eccentricities. Some 
of the largest particles grow from accretion of small bodies, but most remain at their original masses 
(Fig. 3, middle panels). Collisional damping continues until all of the debris-producing particles 
have small eccentricity at t = 4 Myr. Despite their large eccentricities, the largest particles begin 
to grow rapidly (Fig. 3, lower panels). Dynamical friction then reduces their velocities over the 
next 10 Myr. 

To compare these and other results as a function of the starting conditions, we derive the 
damping time from least-squares fits to the luminosity evolution of each model. We use the 
maximum dust luminosity L^ax for these fits. Unless Poynting-Robertson drag is efficient, the 
luminosity from small dust grains is larger than the luminosity contributed by larger objects in the 
main grid. The decline of the dust luminosity is therefore a convenient and relevant measure of the 
damping time. In most models, the dust luminosity declines by a factor of ~ 2 at t = and by a 
factor of ~ 20 at t = lOta- 

Figure 4 shows the variation of the damping time and maximum dust luminosity with x for 
models at 35, 70, and 140 AU. Models with the lowest surface densities are the least luminous 
and have the longest damping times. In systems with x < 10~^, the damping time can exceed the 
stellar lifetime of ~ 1 Gyr. For a given surface density, calculations with larger r^ax have fewer 
small particles. These models therefore are less luminous and have longer damping times. Because 
the collision rate declines with distance from the central star, the luminosity also declines with 
increasing amid- These results indicate that particle collisions in rings with x > 10~^ can produce 
dust luminosities similar to those observed in a real debris disk system with a ring geometry. 
Collisional damping reduces these luminosities by an order of magnitude on timescales of 100 Myr 
or less. 

To test the sensitivity of and to other initial conditions, we vary the initial eccentricity 
Co, the minimum and maximum particle size in the grid {rmin and r„ia-x)'^ the slope 7 of the initial 
mass distribution, Nc{m) oc m~'^; the fragmentation parameters Qc and Sq; and the mass density 
of the particles for calculations at Umid = 35, 70, and 140 AU. Because the dust luminosity is 
proportional to the collision rate, models with larger numbers of small particles are more luminous 
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than models with fewer small particles. These models also damp more rapidly. Thus, models with 

7 > 1 and r^in < 10 cm are more luminous and damp more rapidly than models with 7 < 1 and 
i^min > 10 cm. Models with 'weaker' bodies - smaller Qc and - or with larger, less dense bodies 
are also more luminous and damp more rapidly. These models yield a simple relation between the 
damping time and the luminosity of bodies in the grid at t = 0: 

log {td ■ Ld,o/L^) ^ 0.4 + 1.33 log (a^id/35 AU) . (2) 

The maximum dust luminosity at f = is 

log {td ■ Lmaxfi/L^) ~ 2.1 + 1.33 log (a„jd/35 AU) . (3) 

The scatter in the coefficients is it 0.1 for a wide range of initial conditions. At late times, the dust 
luminosity follows a power law, 

log Lraa./L. OC (t/trf)-^-^^^'^-"^ . (4) 

For the luminosity evolution of the particles in the main grid, we derive 

log Ld/L, OC (iA<i)-'-°^°-' • (5) 

We identify only two exceptions to these relations. In calculations with large eo or weak bodies, 
collisions can remove most of the initial mass on timescales shorter than the relevant damping time. 
Figure 5 shows the ratio of final mass in the ring Mj to the initial mass Mj as a function of cq for 
planetesimals with Sq = 2 x 10^ erg g~^. For modest perturbations with eo < 0.03 (ai/35 AU)^/^, 
collisions damp particle velocities before the ring loses a significant amount of mass to small particles 
with Ti < Tmin- Thcsc modcls follow the Ld — td relations. For strong perturbations with eo > 0.04- 
0.06 {ai/35 AU)^/^, erosion removes a significant amount of mass from the ring before damping 
reduces particle velocities. These models have 25%-50% longer damping times; the asymptotic 
decline of the luminosity is identical to models with smaller perturbations, Ld oc for the large 
particles and L^ax for the small particles. 

In calculations with very large bodies (rmax ^ 500 km), viscous stirring and dynamical friction 
by the largest bodies in the grid can counteract collisional damping. The dust luminosity can then 
remain large for timescales of ~ 100 Myr or longer. Because collisional cascades in rings with 
embedded planets - and no stellar flyby - behave in a similar fashion, we plan to consider these 
models in more detail in a separate paper. 

3.2. Cascades with Gas and Poynting-Robertson Drag 

Gas drag can be an important factor in planetesimal evolution (e.g., Nakagawa et al. 1983; 
Wetherill &; Stewart 1993). Interactions between solid bodies and the gas damp particle velocities 
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and remove particles from the model grid . Both processes scale with the local gas density, pg^i. To 
understand how these processes affect the dust luminosity and the damping time, we calculated a set 
of ring models for different initial values of the gas-to-dust mass ratio, Xg = Tig/Tid- The interstellar 
medium has Xg = 100; observed limits in debris disk systems are uncertain (Zuckerman et al. 1995; 
Dent et al. 1995; Lecavelier des Etangs et al. 1998; Greaves, Coulson, k, Holland 2000a; Thi 
et al. 2001). The appendix summarizes our treatment of particle removal and velocity damping. 
We held the gas density fixed with time; Xg therefore increases with time as fragmentation and 
gas drag remove solid bodies from the grid. Calculations with Xg = 1 are indistinguishable from 
models without gas. Gas drag modifies the evolution of solid bodies for Xg > 10. The following 
paragraphs describe the evolution for Xg = 10 and Xg = 100. 

Figure 6 shows the mass and velocity distributions at several times during the evolution of a 
model with Xg = 10. The initial conditions are identical to the model of Figure 1: small bodies 
with rmin = 10 cm to 10 m have the same initial eccentricity, cq = 0.02, and a total mass of 0.67 
Mq, in 16 annuli at amid = 35 AU from a central 3 Mq star. Interactions between solid particles 
and the gas have a modest impact on the evolution of this model. After 1 Myr, gas drag has swept 
less than 0.01% of the initial mass through the innermost annulus and out of the grid. Compared 
to a model without gas drag, orbital eccentricities are ~ 5% to 10% smaller and the largest bodies 
are ~ 10% to 20% larger (Figure 6, middle panels). The amount of mass lost to small particles is 
~ 24%, compared to ~ 25% in a model without gas drag. After 3 Myr, gas drag is responsible for 
only ~ 0.02% of the total mass loss. Orbital eccentricities are now ~ 20% smaller, and the largest 
bodies are ~ 40% larger (Figure 6, lower panels). This model reaches runaway growth earlier than 
a model without gas drag; the final mass distributions of models with and without gas drag are 
indistinguishable (Kenyon &; Luu 1999). 

Interactions between solid bodies and the gas are much more important when Xg = 100 (Figure 
7). After 0.5 Myr, gas drag reduces orbital eccentricities of the smallest bodies by a factor of four, 
and the largest bodies begin to grow by mergers (Figure 7, middle panels). The mass lost from 
gas drag is small, ~ 0.1% of the initial mass, compared to the mass lost from erosive collisions, 
~ 22.5%. By 1 Myr, gas drag, dynamical friction, and collisional damping cool the small bodies 
well before the largest objects grow to sizes of 100 m (Figure 7, lower panels). The small orbital 
eccentricities of this model yield an early runaway growth phase and the production of several 
Pluto-sized objects on timescales of 10-30 Myr (Kenyon & Luu 1999). The mass of this model at 
1-3 Myr is ~ 3% larger than models without gas drag at similar phases; this difference does not 
yield significantly larger objects after 10-30 Myr compared to models without gas drag. 

Despite their more rapid evolution, models with gas drag follow nearly the same luminosity- 
damping time relation as models without gas drag. We derive log (ta ■ Lmax,o) = C -|- 1.33 log 
(ami(i/35 AU) for these models. Damping is more efficient when gas drag is important; larger gas- 
to-dust ratios imply smaller values of C. Thus, the damping times derived from models without 
gas drag are upper limits to the lifetime of a collisional cascade produced by a stellar fly-by. We 
also derive the same power-law decline of dust luminosity with time, log Lmax/L-t, oc (i/id)~^'^^) as 



- 11 - 



in models without gas drag. 

To complete this portion of our study, we consider mass loss from Poynting-Robertson drag 
for various values of the stellar luminosity and the particle mass density. Mass loss from Poynting- 
Robertson drag is less than 1% of losses from fragmentation and gas drag unless the stellar lu- 
minosity exceeds 100 Lq (p^/l-S g cm~^). Damping of the horizontal velocity is also negligible. 
For A- type stars with < 100 Lq, Poynting-Robertson drag does not modify our conclusions 
unless the mass density of the mid-sized bodies is much smaller than 1 g cm~^. Calculations with 
low density bodies yield smaller damping times than implied by equations (1) and (2). Although 
Poynting-Robertson drag generally is not important for the large grains considered here, removal 
of smaller grains can modify the optical depth (see the appendix). We plan to investigate this 
possibility in a future study. 



3.3. Cascades in Disks 

The ring calculations in §3.1 and §3.2 yield several interesting conclusions about collisional 
cascades in planetesimal disks. Stellar flybys can produce long-lived cascades in low mass disks, but 
the dust luminosity produced in these cascades is small, Lmax/Li, ^ 10~^ for x ^ 10~^. Cascades in 
rings with 1% or more of the minimum mass solar nebula produce larger dust luminosities but have 
short observable lifetimes. The damping times in these models are ~ 1—10 Myr. These timescales 
are short compared to the ages of the oldest debris disk systems, ~ 100-500 Myr (Habing et al. 
1999; Song et al. 2000a; Spangler et al. 2001). 

To test these conclusions in more detail, we consider more complete disk models. The disk 
consists of 64 annuli with Aaj/cj = 0.025 and extends from aj„ = 30 AU to 

^out — 150 AU. As in 

§3.1, the disk contains bodies with rmin = 10 cm to rmax and a mass density of 1.5 g cm~^. All of 
the bodies have the same horizontal velocity with respect to a circular orbit; the initial eccentricity 

1/2 

in each annulus is Cj = eia- where ei is the eccentricity in the innermost annulus. This condition 
yields a constant initial horizontal velocity in the disk. The mass spacing factor is Si = mfe+i/mfe 
= 2 for all batches; the initial mass distribution has equal mass per batch. We adopt the Wetherill 
& Stewart (1993) fragmentation algorithm with Qc = 5 x 10'' erg g~^, Sq = 2 x 10^ erg g~^, and 
= 1 m s-^ 

Fi gures 8^12 outline the evolution of a model with x — 0.1, r^aa; — 

10 m, and ei = 0.02. The 

initial surface density is Sj = 3 g cm" -2 (oi/l AU)-3/2; the initial mass in solid material is 10 M^. 
At the inner edge of the disk, collisions rapidly damp orbital eccentricities on timescales of 1 Myr 
or less. The largest bodies double their radius in the first 1 Myr and quadruple their radius in the 
first 10 Myr. Because collision cross-sections scale with disk radius as ^ oc a~^, collisional damping 
proceeds very slowly at the outer edge of the disk. Growth is slow. Collisional damping starts 
to affect mid-sized bodies in the outer disk at 1 Myr, when bodies at Oj < 40 AU are completely 
damped and starting to grow rapidly (Figure 8; middle panels). Gravitational focusing leads to 



- 12 - 



faster growth once mergers dominate erosive collisions. Thus, large bodies at 30 AU have entered 
the runaway growth phase when 10 cm bodies at 150 AU are just starting to grow (Figure 8; lower 
panels) . 

Figure 9 shows the sensitivity of collisional damping and of the growth rate of the largest 
objects to position in the disk. The initial eccentricities of the smallest particles follow a power- 

1/2 3/2 

law, oc ; the scale height of these bodies above the disk midplane is thus H (x . The 
eccentricity of the smallest particles in the grid damps by a factor of ~ 2 at = 30 AU in 1 Myr, 
when small particles at the outer edge of the grid have almost all of their initial orbital eccentricity. 
At 10 Myr, damping has reduced the orbital eccentricities of all small particles by factors of 2 or 
more. The eccentricities of the smallest particles then vary with disk radius as cx a°'^^. The 
power-law slope of this relation falls to 0.8 at 50 Myr. This behavior yields a very steep variation 
of the scale-height with radius, H oc a\'^. 

The merger rate is more sensitive to disk radius. For t < 10 Myr, growth rates are linear, 
because gravitational focusing factors are small. The radius of the largest body is a smooth function 
of disk radius, Vi^rnax ^ ol^^ with 7 > 1.5. Once particles in the inner disk reach sizes of ~ 1 km, 
gravitational focusing factors increase and growth becomes non-linear. Annuli in the runaway 
growth phase depart from a smooth power-law variation of n^rnax with Oj, as indicated by the sharp 
rise in Vi rnax 

for Qi < 40 AU at 50 Myr. 

The total dust luminosity of this model follows the evolution derived for the ring models in 
§3.1 (Figure 10). At t = 0, the maximum dust luminosity of Lmax/Li, ~ 6 x 10~^ is comparable 
to luminosities observed in [3 Pic and other luminous debris disk systems (Augereau et al. 1999a; 
Barrado y Navascues et al. 1999; Habing et al. 1999; Song et al. 2000b). The dust luminosity 
decreases with time and reaches Lmax/Li, ~ 2 x 10~^ at t = 50 Myr. This luminosity is close 
to luminosities observed in a Lyr and other faint debris disk systems (Kalas 1998; Spangler et 
al. 2001). Equation (1) describes the time-evolution of the luminosity remarkably well. The 
luminosity-damping time relation is 

log td ■ L^ax,o/L^ = 3.16 ± 0.05. (6) 

At late times, 

Lmax/L. « 4 X 10-3 (Vtrf)-i-°^±°-°' . (7) 
Kalas (1998) derives L oc for known debris disk systems with reliable ages. 

To visualize the time evolution of the luminosity for this disk model, we assume a pole-on 
system where each annulus in the grid has a dust opacity derived in the geometric optics limit from 
the number density (§A.6). The luminosity per unit surface area radiated by each annulus is then a 
simple function of the opacity and the dust scale height (Kenyon et al. 1999). We adopt the scale 
height for the smallest particles in the grid. This procedure yields the bolometric surface brightness. 
If the particle albedo is large and independent of grain size, the predicted surface brightness should 
be close to the optical or near-IR surface brightness of scattered light. 
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Figure 11 shows how the radial surface brightness profile changes with time. At the start of our 

calculations, models with a power-law surface density distribution Ej produce a power-law surface 

3/2 7 /2 

brightness distribution For Sj oc a- , we derive Ii cc a- . For t < \ Myr, collisions damp 
particle eccentricities in the inner disk, ri ~ 30-60 AU. Smaller collision velocities produce less dust; 
the surface brightness fades. At f = 1 Myr, the surface brightness is nearly constant for rj = 30-60 
AU and gradually approaches the initial ij oc at > 100 AU. As the calculation proceeds, 

collisions damp particle velocities in the outer disk. After i ~ 10 Myr, the surface brightness in the 
inner disk declines by nearly two orders of magnitude; the outer disk fades by ~ 30%. All of the 
disk then fades by roughly an order of magnitude over the next 20-40 Myr. When we terminate 
the calculation at 50 Myr, the surface brightness distribution is much more shallow than the initial 

j^/2 7/2 

distribution, with li oc a- instead oi li oc a- 

Figure 11 illustrates how coagulation concentrates mass into large objects with small optical 
depth. At 50 Myr, the radial gradient in the surface density of solid material is similar to the 
initial gradient, with S oc a~^'^^. Because material in the inner disk has more frequent collisions 
than material in the outer disk, the inner disk has more material in large objects with r,j > 100 
m than the outer disk. The surface density of small objects thus grows with radius; for t > 10 
Myr, 

'^small flj' for < 10 cm. There is some evidence for this behavior in real systems: the 
particle disks in BD-h31°643 (Kalas & Jewitt 1997) and HD 141569 (Weinberger et al. 1999) 
have regions where the number density of small particles appears to increase in radius. We plan 
addititonal calculations to see if this feature of the calculations is characteristic of flyby models or 
all planet-forming particle disks. 

Figure 12 shows disk images at six times in the evolutionary sequence. An animation of the 
disk evolution is available in the electronic version of this paper. To construct predicted images from 
the surface brightness profiles, we assume a simple power law relation to convert surface brightness 
to pixel intensity. We simulate an observation by adding counting noise to the intensity of each 
pixel, but we do not convolve the image with a 'seeing disk.' These approximations ignore radiative 
transfer effects in the disk and the variation of disk flux with wavelength, but the resulting images 
accurately illustrate the physical behavior of the dust luminosity with time. 

To investigate the sensitivity of these results to the starting conditions, we consider several 
models with different values for ei and x and several models with gas drag. Full disk models require 
more computer time than ring models; this study is therefore more limited than the parameter study 
of §3.1. As long as erosive collisions do not exhaust the supply of mid-sized bodies before collisional 
damping can reduce particle eccentricities, models without gas and Poynting-Robertson drag yield 
the luminosity-damping time relation in equation (6) to an accuracy of 10% or better independent 
of the initial conditions. Low mass disks are less luminous and longer lived than massive disks. 
Calculations with gas drag {Xg = 100) evolve ~ 30% faster than models without gas drag. Equation 
(5) thus provides an upper limit to the damping time for a disk with initial luminosity L^axfi/Li,. 

The evolution of the scale height is the main difference between models with and without gas 
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drag. When the gas density is low {X^ < 1-3), gas drag does not reduce significantly the orbital 
eccentricities of the smallest particles in the grid. Collisional damping and dynamical friction 
dominate the velocity evolution; the scale height of small particles exceeds the scale height of large 
particles. When the gas density is larger (X^ > 10), the small particles are more closely coupled to 
the gas and have more circular orbits. The mid-sized particles then have the largest scale heights. 
If the gas and small particles are well-mixed perpendicular to the disk midplane, the dust has a 
large scale height and a large dust luminosity (Kenyon Sz Hartmann 1987). If mixing between 
dust and gas is poor, small dust particles settle to the midplane; the dust luminosity is then small. 
Detailed hydrodynamical calculations are necessary to derive the outcome of a collisional cascade 
in a planetesimal disk with large amounts of gas. 

Full disk calculations yield the same behavior with initial eccentricity as ring models. For 
annuli with eg < 0.05 (a^/SS AU)^/^, collisional damping reduces particle velocities before the disk 
loses more than half of its initial mass. The disk loses more than 95% of its initial mass for strong 
perturbations with eo > 0.05 (ai/35 AU)^/^. In these strong perturbations, collisional damping 
reduces particle velocities sufficiently to allow growth by mergers and the formation of 1 km or 
larger bodies. The timescale for planet formation is then very long, ~ 1 Gyr or longer. 

Our results indicate that complete disk models evolve on faster timescales than ring models (see 
also Kenyon & Bromley 2001). In ring models, planetesimals at the edges of the grid interact with 
fewer annuli than planetesimals in the middle of the grid. Collisional damping is less effective with 
fewer collisions. Partial disk models thus damp more slowly and lose more mass than complete disk 
models. At 35 AU, complete disk models lose ~ 20% less mass than ring models. This difference 
rises to ~ 50% at 140 AU. More effective damping and less mass loss yield larger bodies on shorter 
timescales in complete disk models. 

3.4. Limitations of the Models 

Statistical simulations of planetesimal growth in a dusty disk have a long history, with well- 
known limitations and uncertainties (Wetherill 1980; Greenberg et al. 1984; Davis et al. 1985; 
Barge & Pellat 1991; Spaute et al. 1991; Lissauer & Stewart 1993; Wetherih & Stewart 1993; 
Stern &Colweh 1997a; Weidenschilling et al. 1997; Kenyon & Luu 1998; Inaba et al. 2001). The 
particle-in-a-box formalism assumes a homogeneous spatial distribution of solid bodies with small 
velocities relative to the local circular velocity. This assumption is good until runaway growth 
produces a few massive bodies not distributed uniformly in space (Kokubo & Ida 1996). Our 
simulations never reach this limit. Our calculations also require a specific algorithm for collision 
outcomes (§2 and the appendix). As long as particle collisions produce damping and some debris, 
our derived relations for the evolution of luminosity as a function of initial disk mass are remarkably 
independent of the details of the collision algorithm. The collisional damping prescription is more 
important; changing the coefficients in the algorithm (equations AS and A9 of the appendix) changes 
the damping times by a similar factor. We have not investigated the sensitivity of the results to 
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the form of the damping equations. We have tested our prescription for velocity evolution in the 
low velocity limit, where the relative collision velocity is small compared to the Hill velocity (Barge 
k Pellat 1990; Wethcrill & Stewart 1993; Kenyon & Luu 1998; Stewart & Ida 2000; Kenyon 
& Bromley 2001). The statistical approach is invalid in this limit. Ida (1990), Barge & Pellat 
(1990), Ida & Makino (1992, 1993), and Wetherih & Stewart (1993) have described solutions to 
this problem; we use the Ida & Makino (1992) approach. Tests with other approaches yield the 
same results. 

The implementation of the algorithms produces other limitations. Our multiannulus calcula- 
tions eliminate many of the constraints associated with calculations of a single accumulation zone 
(Wetherih 1990b; Wetherih & Stewart 1993; Weidenschilling et al. 1997; Kenyon & Luu 1999). 
We can treat collisions and gravitational perturbations between particles in adjacent annuli, and 
the drift of particles through adjacent annuli due to gas drag and Poynting-Robertson drag. The 
current calculations do not include algorithms for orbital migration of planetesimals or the evo- 
lution of gas in the disk (including gas accretion), but these processes are not important for the 
situations we investigate here. 

Uncertainties resulting from the finite mass and spatial resolution should be small. Calculations 
with coarse mass resolution within an annulus lag higher resolution simulations by 10% to 20% but 

yield nearly the same mass and velocity distributions (see Kenyon & Luu 1998, 1999). Calculations 
with finer spatial resolution are indistinguishable from the simulations described above at the 5% 
level. These differences grow with the mass of the largest object in the grid but are not important 
for the relatively small mass objects we consider here. 

To compute the dust luminosity, we assume two limits to the dust opacity of small objects. 
We make these assumptions to avoid computing the time evolution of small particles, which have a 
negligible impact on the evolution of larger planetesimals in the grid. In most planetesimal disks, 
the large opacity limit is appropriate, because the timescales for Poynting-Robertson drag are ~ 1 
Myr or longer at 30 AU and ~ 25 Myr or longer at 150 AU. These timescales are long compared 
to the damping times of massive, luminous planetesimal disks. For more luminous central stars 
with Li, ~ 100 Lq, the small opacity limit is appropriate, because Poynting-Robertson drag can 
remove small grains on timescales much shorter than 1 Myr at 30 AU. We plan to address these 
uncertainties in more detail in a future paper. 

In our complete disk models, we adopt a shallow variation of the initial eccentricity of plan- 

1/2 

etesimals with disk radius, cq oc a/ . The ratio of the initial horizontal velocity to the shattering 
velocity of a dust grain is then independent of Oj, which allows us to investigate collisional damping 
without worrying about the importance of different fragmentation algorithms. This approximation 

is reasonable for moderately strong perturbations where the passing star has a periastron of ~ 600 

5 /2 

AU. Closer stellar fiybys produce stronger perturbations with cq (x a- (Ida et al. 2000; Kalas 
et al. 2000; Kobayashi &; Ida 2001). In a strong perturbation, planetesimals in the outer disk 
are depleted more rapidly than those in the inner disk. Even if the outer disk loses nearly all of 
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its initial mass, our calculations indicate that collisional damping should reduce particle velocities 
and allow planet formation on long timescales. We plan to investigate this possibility in a future 
paper. 

The speed of the computer is a final limitation to these calculations. Although larger test 
simulations are possible, practical considerations limit calculations of model grids to 64-128 annuli. 
For reasonable spatial resolution, the maximum extent of our grid is a factor of 5-10 in distance 
from the central star. We are thus 1-2 orders of magnitude from constructing model grids of 
complete solar systems. For a given set of algorithms, the timescales and dust luminosities derived 
here should be within 10% to 50% of those derived from more extensive models. Limitations on the 
extent of the radial grid are more important for addressing the interfaces between (i) gas giants and 
terrestrial planets and (ii) gas giants and the Kuiper Belt in our solar system. Similar problems 
occur in evolutionary calculations of other types of circumstellar disks. Faster computers should 
resolve these difficulties in the next few years. 

4. DISCUSSION AND SUMMARY 

Stellar flybys can plausibly produce collisional cascades in a planetesimal disk. When the mass 
of solid material in the disk is 1% or more of the minimum mass solar nebula (~ 1-100 of solid 
material at 30-150 AU), flybys which increase planetesimal velocities close to the shattering limit 
can produce dust luminosities similar to those observed in known debris disk systems. The dust 
luminosity correlates with the disk mass. Collisional damping in these disks is very efficient; dust 
production decreases dramatically on timescales of 1 Myr or less. The most luminous and most 
massive disks damp the fastest. These conclusions are independent of many uncertainties in the 
calculations, including the initial mass distribution, the fragmentation parameters, and the mass 
density of planetesimals. 

Our calculations indicate two outcomes for a collisional cascade induced by a moderately close 
stellar fiyby. If the fiyby produces relative velocities which exceed the strength of a planetesimal, 
disruptive collisions nearly exhaust the supply of planetesimals before collisional damping becomes 
effective. Most planetesimals are eventually ground to dust and planet formation is difficult (see 
also Kobayashi Sz Ida 2001). For icy planetesimals with 5*0 = 2 x 10® erg g~^, this limit is cq > 
0.05 (ai/35 AU)V2 

When a flyby produces a smaller perturbation in the disk, collisional damping reduces plan- 
etesimal velocities on relatively short timescales. This process allows the disk to retain a significant 
mass in large planetesimals, which collide and merge to form planets. Our ability to detect the 
collisional cascade observationally depends on the disk mass. Collisional cascades in disks with 
~ 1% or more of the minimum mass solar nebula produce dust luminosities, L„iax/Li, > 10~^, 
comparable to those observed in debris disk systems (Kalas 1998; Spangler et al. 2001). 

Collisional damping is a new feature in models of the evolution of collisional cascades in 



-17- 



planetesimal disks. Previous models assumed that collisions preserve the eccentricity and inclination 

distributions of planetesimals in the disk (e.g., Mouillet et al. 1997; Augereau et al. 2001, see 
also Artymowicz 1997, Lagrange et al. 2000, and references therein). Dust grains produced by 
collisions thus had the same orbits as their progenitors. In our models, 'collisionless' systems 
with long collision timescales yield optical depths and dust luminosities much smaller than those 
observed (Figures 4 and 10). Systems with larger optical depths and dust luminosities have short 
collision times and short damping timescales. This conclusion is an extension of earlier work on 
the debris disk system HR 4796A and the Kuiper Belt in our solar system, where disks with masses 
of solid material comparable to or larger than the minimum mass solar nebula allow the growth 
of Pluto-sized objects, which stir the velocities of smaller objects to the shattering velocity and 
produce observable collisional cascades (Kenyon & Luu 1999; Kenyon et al. 1999, see also Stern 
& Colwell 1997a,b). 

Following a flyby, equation (1) fits the evolution of the dust luminosity to remarkably high 
accuracy. Because damping and dust production depend on the collision rate, collision physics 
sets the form of this relation. To derive the limiting form of equation (1), we assume an ensemble 
of identical planetesimals. The coagulation equation is dn/dt oc (§A.2). The dust production 
rate is dm/dt oc m?'n? oc v? where m is the mass of a planetesimal. Thus, the disk mass and the 
number density of planetesimals fall with time, oc t~^. Because the disk has constant surface 
area, the surface density follows the same relation, S oc The left panels of Figure 3 verify this 
relation for the damping time, td oc S^^. The luminosity of the planetesimals is proportional to 
the disk mass, which yields L oc t^^; our models follow this relation faithfully (Figure 2 and 10). 
The dust luminosity depends on the mass distribution of planetesimals for < Vmin, where r„jj„ 
is the minimum planetesimal radius in the grid. The appendix describes two limits on this mass 
distribution. The maximum dust mass assumes that Poynting-Robertson and other drag forces 
do not remove particles from the grid; the minimum dust mass assumes that drag forces remove 
particles efficiently from the grid. In both limits, the dust mass is proportional to the disk mass 
(Figure 3, right panels). The dust luminosity thus declines with time, 

Lmax,min OC t . (8) 

Kalas (1998) derives the same relation from observations of debris disks. Our detailed models for 
complete disks yield this relation; because planetesimals have fewer annuli to interact with, the 
relation for ring models is steeper. 

Spangler et al. (2001) derive a different relation for the time-variation of luminosity in a 
collisional cascade, oc t~^. They assume that the dust clearing time is shorter than the age of the 
system and conclude that the dust luminosity is proportional to the instantaneous dust production 
rate, dn/dt. However, the dust luminosity is proportional to dn/dt only if planetesimal collisions 
produce dust which escapes the system on timescales shorter than the collision time. Because 
Poynting-Robertson drag and gas drag are slow processes in the outer disk, radiation pressure is 
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the only mechanism capable of removing particles from a large debris disk on short timescales^. 

Particles with sizes of 10 fiui or larger are unaffected by radiation pressure. Planetesimal collisions 
thus build up a reservoir of dust grains with sizes of 10 ^m or larger, while radiation pressure 
ejects smaller grains (Burns, Lamy, & Soter 1979; Backman &; Paresce 1993; Artymowicz 1997). 
For the power law cumulative mass distribution appropriate for collision debris, Nc{m) oc rn~^'^^ 
(Dohnanyi 1969; Davis et al. 1985; Wetherill & Stewart 1993; Kenyon & Luu 1999), the optical 
depth of grains ejected by radiation pressure is a factor of 10 or more smaller than the optical 
depth of larger grains. The larger grains in an annulus thus make a larger contribution to the 
dust luminosity than smaller grains accelerated away from the annulus. The decline of the dust 
luminosity then depends on the rate of decline for the total disk mass, Lmax oc t~^. 

Our results suggest modifications of the picture for the formation and early evolution of the 
Kuiper Belt. Ida et al. (2000) propose a stellar flyby to explain the structure of the outer Kuiper 
Belt in our solar system. In their model, a stellar encounter with a distance of closest approach of 
~ 100 AU excites the large orbital eccentricities and inclinations observed for 50-500 km objects at 
distances of 30-50 AU from the Sun. Collisions between these objects lead to a collisional cascade 
which erodes the mass of the inner Kuiper Belt over time. In our calculations, collisional damping of 
1-100 m planetesimals accompanies mass erosion; dynamical friction then damps the eccentricities 
and inclinations of the largest objects on short timescales, ~ 1-10 Myr for a minimum mass solar 
nebula model and ~ 100 Myr for a disk with a mass of solid material equal to the mass of the 
current Kuiper Belt. Our models predict damping times of 5 Gyr or longer only when the disk has 
a small mass in solid objects with radii of 1 km or less. Coagulation models which can produce 
50-500 km Kuiper Belt objects at 30-40 AU leave most of the initial mass in smaller objects with 
radii of 0.1-10 km (Kenyon 2002; Kenyon & Luu 1999, see also Stern & Colweh 1997a,b). Thus, if 
a stellar encounter is responsible for the large orbital eccentricities observed in Kuiper Belt objects, 
our models suggest that this encounter occurred after some other process depleted the Kuiper Belt 
of 0.1-1 km objects. We plan to investigate the idea that stirring by one or more Pluto-sized objects 
embedded in the Kuiper Belt can explain this depletion. 

Kalas &L Jewitt (1995, see also Kalas et al. 2000; Kalas, Deltorn, Sz Larwood 2001; Larwood & 
Kalas 2001) propose a stellar flyby model to explain structures observed in the (3 Pic debris disk. In 
their picture, a recent, ~ 10^ yr, stellar encounter produces an asymmetric ring system and other 
dynamical features in the disk (Kalas et al. 2000). The highly perturbed orbits of planetesimals in 
the disk lead to a collisional cascade and copious dust production. Although we cannot calculate 
planetesimal evolution in the ~ 1000 AU P Pic disk with adequate spatial resolution, our 30-150 AU 
disk model in §3.3 provides a reasonable first approximation to the evolution of a collisional cascade 
induced by a moderately close stellar flyby in the inner portion of a large disk. The maximum dust 



^Poynting-Robertson and gas drag can efficiently remove small particles from a planetesimal disk at distances of 
1-5 AU from the central star. Relatively rapid dust removal might account for the large inner 'holes' observed in 
some debris disk systems (see also Backman & Paresce 1993; Artymowicz 1997; Koerner et al. 1998; Lagrange et 
al. 2000). 
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luminosity of Lmaxfi/Li, ~ 4 x 10~^ is close to the bolometric luminosity ratio L^/Li, 3 x 10~^ 
observed in /3 Pic (Backman & Paresce 1993). The damping time of the model, ~ 3 x 10^ yr, 
indicates that collisions do not wash out either the observed ring structure or the radial variation 
of eccentricity and inclination at distances of 100-1000 AU on timescales similar to the apparent 
dynamical age of ~ 10^ yr. However, our models produce significant coUisional damping at SO- 
SO AU on short timescales. Damping produces an apparent hole in the disk, with smaller dust 
production than annuli in the outer disk (Figures 11-12). In principle, our models can produce 
changes in the slope of the radial surface brightness distribution similar to that observed in the 
(3 Pic disk (Golimowski, Durrance, & Clampin 1993; Augereau et al. 2001). Radiative transfer 
calculations similar to those of Kenyon et al. (1999) will allow us to make a detailed comparison 
between our models and observations. 

Associating stellar flybys with other debris disk systems seems less practical. The damping 
time-luminosity relation (equation (6)) implies that all collisional cascades induced by stellar flybys 
have luminosities which decline below L^ax/Li, ~ 10~^ on timescales of 100 Myr or less. Source 
statistics suggests that nearly all stars are born with a disk (Lada 1999); many main sequence 
stars also have disks (Habing et al. 2001). Several known debris disk systems have ages of 500 
Myr or more (Song et al. 2000b; Spangler et al. 2001). Recent flybys for all objects are unlikely; 
the chance probability of a stellar encounter which passes within ~ 500 AU of a nearby field star 
is 1% or less in ~ 100 Myr (Frogel & Gould 1998; Garcia-Sanchez et al. 1999). We suspect that 
stirring by embedded planets or a binary companion induces collisional cascades in most debris 
disk systems. We plan future papers to address this issue. 

Our results suggest caution when interpreting observed relations between dust luminosity and 
age in debris disk systems. Most evolutionary processes in a planetesimal disk scale with the orbital 
period. The factor of 4-6 range in the stellar mass among known debris disk systems implies a 
range of 2-2.5 in evolutionary timescales for systems with the same initial disk mass and disk radius. 
In an ensemble of stars with a range of temperatures and luminosities, Poynting-Robertson drag 
will cause lower dust luminosities in disks with luminous stars compared to less luminous stars. 
Poynting-Robertson and gas drag may also reduce disk luminosities below current detection limits 
on shorter timescales for luminous stars than for less luminous stars. This behavior suggests that 
the timescalc for disk evolution should correlate with the nature of the central star. Song et al. 
(2000b) comment on this effect in their data for 10-15 debris disk systems. Although the trend 
described by Song et al. (2000b) is encouraging, current samples are too small and the selection 
effects too uncertain to make robust statements about timescales for disk evolution. Larger samples 
with smaller detection limits may provide enough data to test trends predicted by the models in 
detail. 

Finally, multiannulus accretion codes are an important step towards understanding the for- 
mation and evolution of planetary systems (see also Spaute et al. 1991; Weidenschilling et al. 
1997). Our calculations demonstrate that 'complete' disk models with 64 or more annuli evolve on 
faster timescales than ring models with 16 or fewer annuli. Planetesimals at the edge of a model 
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grid interact with fewer annuli than planetesimals in the middle of a grid. Colhsional damping, 
dynamical friction, and runaway growth thus proceed more rapidly in a complete disk model than 
in a ring model. Calculations of fragmentation, gas drag, and Poynting-Robertson drag arc also 
more accurate in a multiannulus grid. This difference is important for constructing accurate models 
for the formation of Jupiter and the Kuiper Belt, where observational constraints limit the mass 
and timescale available for planet formation (e.g., Lissauer 1987; Bailey 1994; Pollack et al. 1996; 
Kenyon & Luu 1999; Bryden, Lin, & Ida 2000) 
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ics Origin-2000 'Alhena' through funding from the NASA Offices of Mission to Planet Earth, Aero- 
nautics, and Space Science. Advice and comments from M. Geller, M. Kuchner, and an anonymous 
referee greatly improved our presentation. 
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A. APPENDIX 
A.l. Overview 

We assume that planctesimals are a statistical ensemble of masses in N concentric, cylindrical 
annuli with width Aoj and height Hi centered at radii from a star with mass and luminosity 
L^. The particles have horizontal hik{t) and vertical Vik{t) velocity dispersions relative to an orbit 
with mean Keplerian velocity Vxi (Lissauer &: Stewart 1993; Stewart k. Ida 2000). We approximate 
the continuous distribution of particle masses with discrete batches having an integral number of 
particles nik{t) and total masses Mik{t). The average mass of each of M mass batches, mik{t) = 
Mik{t)/nk{t), evolves with time as collisions add and remove bodies from the batch (Wetherill & 
Stewart 1993). 

Kenyon & Luu (1998, 1999) describe our approach for solving the evolution of particle numbers 
and velocities for a mass batch A; in a single annulus i during a time step St. Kenyon & Bromley 
(2001) describe our solution for the velocity evolution from elastic collisions for a set of annuli. 
Here we describe collision rates and velocity evolution from inelastic collisions for a set of annuli. 
We also include our treatment of particle loss and velocity evolution for gas drag and the Poynting- 
Robertson effect. 

A. 2. Coagulation Equation 

The coagulation equations for a particle in mass batch k of annulus i colliding with another 
particle in mass batch I of annulus j are 

Sriiik' = St [eijkiAijkinikriji — niik'Ai/jk'inji] + Siii'k'j — Sni'k',gd — ^^i' k' ,prd (Al) 



SMiiki = St niiiki [eijkiAijkinikriji — ni/kiAi/jkunji] + Srriiik'j — Srriiiki^gd — Suiiik/^prd (A2) 

where Aijki is the cross-section, Cijki = 1/2 for i = j and k = I, and Cijki = 1 for k ^ I and 
any The terms in A1-A2 represent (i) mergers of and rriji into a body of mass mj/fc/ = 
rriik + rriji — m^ijki, (ii) loss of nii'k' through mergers with other bodies, (iii) addition of nii'k' from 
debris produced by the collisions of other bodies, and (iv) loss of mj/^/ by gas drag and Poynting- 
Robertson drag. Kenyon &; Luu (1999) outlines our calculation of the mass lost to small fragments, 
iTT-e,ijki- The second term in A1-A2 includes the possibility that a collision can produce debris but 
no merger (rebounds; see Davis et al. 1985; Barge & Pellat 1993; Kenyon & Luu 1999, and 
references therein). 
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Following Weidenschilling et al. (1997), particles in adjacent annuli may collide if their orbits 
approach within 2.4 times their mutual Hill radius Rh- The 'overlap region' for these inelastic 
collisions is 



Oijkl,in = 2ARh + 0.5{wik + Wji) - \ai - aj\ ; (A3) 

where is the radial extent of the orbit of particle k with orbital eccentricity in annulus i 
(Kenyon h Bromley 2001): 



Wik={ , „ „ „ /A„ M/4 A„ (A4) 



(Aoj + efeai)(efeai/Aai)^/^ efeOj > Aa 



We calculate the 'gravitational range' of the largest bodies - Rg^ik = KiaRn^ijki + 2aejfe (Wetherill 
&; Stewart 1993) - where Ki = 2\/3 and RH,ijkl = [{nT'ik + "ijO/^^ol^^^ the mutual Hill radius 
(see also Kenyon & Luu 1998, 1999). Isolated bodies are the N largest bodies that satisfy the 
summation, ^^^""^ nikRg^ik > Acj (Wetherill & Stewart 1993). These isolated bodies cannot 
collide with one another but can collide with other lower mass bodies in the annulus and all other 
bodies in other annuli. 

We use the particle-in-a-box technique to calculate collision cross-sections, A^jki- Kenyon & 
Luu (1998) describe this approach for a calculation with a single annulus. The technique is easily 
generalized to a multi-annulus calculation by averaging quantities which depend on the radius of 
the annulus, such as the Keplerian velocity or the Hill radius (Spaute et al. 1991; Weidenschilling 
et al. 1997). The cross-section is 

Aijki = oicoii (777 / \ /A — \] ^ijki ^9,m (rik + Vjif , (A5) 

where acoii is a constant (Wetherill & Stewart 1993; Kenyon & Luu 1998), Hijki is the mutual 
scale height, {aij) and (Aajj) are the average heliocentric distance and width for the two interacting 
annuli, Vijj^i is the relative particle velocity, Fg^ij^i is the gravitational focusing factor, and rj^ and 
rji are the particle radii. We adopt the piecewise analytic approximation of Spaute et al. (1991) 
for the gravitational focusing factor, and the two body coUisional cross-sections of Greenberg et al. 
(1991) in the low velocity limit (see also Greenzweig & Lissauer 1990, 1992). 

To calculate i' and k' for each collision, we establish fixed grids of annuli and masses. The 
radial grid has annuli with either a fixed constant size, Aoj = aj_|^i — flj, or a variable width, 
Aflj = Aao • fli . The mass grid in each annulus has particle numbers A; = 1 to = N^ax with 6k = 
mk+i/mk- We set 

5i k < ko 

61 - e ■ {k - ko) ko <k <ki (A6) 

62 k > ki 
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where Si, S2, and e are constants. Grids with a fixed mass ratio between neighboring mass bins 
have e = 0. Grids with e > allow coarse resolution, S ~ 1.4-2.0, where small objects grow slowly, 
and fine resolution, 5 ~ 1.05-1.4, where large objects grow rapidly. When a collision produces n^^ 
bodies with m^/ in annulus i' , we augment either batch I' when m^' < ^m;/m;/_)_i or batch I' + 1 
when rrik' > ^ymj/mi/^. We add mass to annulus j' when f — j' < 0.5 and to annulus j' + 1 when 
f — j' > 0.5, where 

J ^ iruik + jrriji 
niik rriji 

and j' is the integer part of /. A complete cycle through all mass batches and all anmili produces 
new values for n^/ and M^/ in each annulus i' , which yields new values for the average mass per 
bin, ruk' = Mki/riki. This process conserves mass and provides a good description of coagulation 
when Sk is small (Ohtsuki & Nakagawa 1988; Wetherill 1990a; Ohtsuki et al. 1990; Kenyon & 
Luu 1998, e.g.,). Our algorithm to place merged objects in a particular annulus does not conserve 
angular momentum specifically; test calculations conserve angular momentum to better than 1% 
after 10^ timesteps. 



A. 3. Velocity Evolution 

The time evolution of particle velocities from collisional damping is 



— ^ = X] X] ^ ("^ii^li - rnikhjj. - {niik + mji)h]k) h (A8) 



(it!? . 3 = Nl = lraax ^ 

— ^ = X X r ~ ^ii'^ik - i'^ik + mji)v^k) k (A9) 

j=Q 1=0 Pijkl 

where Cm = acoii fijkl ^ijkl Pg,jl Vijkl Fg,i3ki {rik + rjif (Hornung, Pellat, & Barge 1985; Wetherill 
& Stewart 1993; Ohtsuki 1992, 1999). In the second summation, Ijnax = k when i = j; Ij^ax = 
M when i ^ j (see also Kenyon & Luu 1998, 1999). We add a term, fijki, to treat the overlap 
between adjacent zones; /jj^; = 1 when i = j and fij^i < 1 when i 7^ j. The integrals Ig and /j are 
elliptic integrals described in previous publications (Hornung, Pellat, & Barge 1985; Wetherill &: 
Stewart 1993; Ohtsuki 1992, 1999, Stewart & Ida 2000). 



A.4. Gas Drag 

Gas drag can be an important process during the early evolution of planetesimal disks. The 
viscosity of the gas drags small bodies inwards through the disk and damps the motion of small 
bodies relative to a circular orbit. Adachi et al. (1976) and Weidenschilling (1977b) analyze 
the interactions between gas and planetesimals, and develop useful evolution formulae (see also 
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Kary, Lissauer, &; Greenzweig 1993). We adopt the Adachi et al. (1976) prescription. Material 
drifts through an annulus at a rate 

Sol St 

= 2(0.97eife + OMiik + Vik/VK,i)Ty- — , (AlO) 

VK,i td,ik 

where rji is the gas velocity relative to the local Keplerian velocity, VK,i- The characteristic drift 
time is 



365 / V/Vl AU\ /ICr^^cm^^X 

where Cij = 0.5 is the drag coefficient, pg^i is the gas density, and Tx^i is the orbital period (Adachi 
et al. 1976; Wetherill & Stewart 1993). The removal of bodies from annulus i at each time step 
is then 

Sngd,ik _ Saik 



riik 



(A12) 



To specify the gas density, we adopt an initial gas-to-dust surface density ratio, 'Ego^i/T,do,i, 
and set pgo^i = Sgo,i/2i?go,i) where -ffgo,i is the vertical scale height of the gas. We assume a scale 
height appropriate for a gaseous disk in hydrostatic equilibrium: 

Hgo^i = Ho{ai/ao)'^ai , (A13) 

where Hq is the scale height at ag = 1 AU. In most cases, we adopt Hq = 0.05 and q = 0.125 
(Kenyon &; Hartmann 1987). We also assume that the gas density falls with time: 

Pg,i{t) = P90,i e-'^'' (A14) 

where tg is a constant and is usually 10 Myr to 1 Gyr. 

Our gas drag algorithm removes bodies from annulus i + 1 and places them in annulus i. We 
explicitly conserve kinetic energy: the velocity dispersion of annulus i decreases as lower-velocity 
material is dragged inwards from outer annuli. The grid loses material and kinetic energy at the 
inner boundary and gains material and kinetic energy at the outer boundary. To compute the 

amount of material added to the grid, we assume a power-law density distribution, pg^i oc a~°'. 
Equation (A12) is then 5ngd^ik/nik oc a- where we have ignored the slow variation of e, i, 

and r]i/VK,i with Oj. The number of particles added to the last zone in the grid is: 

^'n'gd,Nk = ( ) Sngd,Nk ■ (A15) 
\aN+iJ 

The addition of kinetic energy follows 

(\ a-l-1.5 
SEgd,Nk , (A16) 
^N+lkJ 
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where Sugd^Nk SEgd,Nk the number and kinetic energy of particles dragged inwards from 
annulus N to annulus N — 1. 

We adopt the Wetherill & Stewart (1989) expression for velocity damping due to gas drag: 

dt ~ 2m,fe Pa^'^9,ik^ik^ y^^^) 

where Co = 0.5 is the drag coefficient and Vgi = {Vik{Vik + Vi))^^^ is the mean relative velocity of 
the gas. 



A.5. Poynting-Robertson Drag 

When the luminosity of the central star is large, the radiation field can drag small particles 
towards the central star and can circularize particle velocities relative to a circular orbit. Radiation 
pressure can also eject very small particles from the disk. Burns, Lamy, &; Soter (1979) have 
derived accurate expressions for these processes. Because we do not treat the evolution of the 
smallest grains, we ignore radiation pressure. For Poynting-Robertson drag, the inward drift of 
material is 

SUik / 2 + 3e4 \ VprQpr f.-.r.. 

where Qpr is the Mie scattering coefficient and 

2.53 X 10" / \ 
Vpr = [-r^] ■ (A19) 

The mass density of an individual grain is pik- The rate particles leave an annulus is 



(A20) 



As with gas drag, we scale the number of particles lost by annulus N to calculate the number of 
particles added to this zone: 

^npr,Nk = ( — ^ ) Sripr^Nk (A21) 
\aN+iJ 



For the kinetic energy: 



SE'^^^,= (-^)' 6Ep,,^k. (A22) 
\aN+ikJ 



In these expressions, 5npr^Nk aiid SEp^^Nk are the number and kinetic energy of particles dragged 
inwards from annulus N to annulus N — 1. 
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We also adopt the Burns et al. (1979) expression for horizontal velocity damping from 
Poynting-Robertson drag, 



5f^ 



a 



= _6.25 "^^'7'^' . (A23) 



Poynting-Robertson drag does not change the vertical component of the velocity. 



A. 6. Dust Luminosity and Optical Depth 

To calculate the dust luminosity, we follow Kenyon et al. (1999) and estimate the bolometric 
luminosity reprocessed by solid objects in the disk as L/L^, = tCI/Att, where r is the radial optical 
depth and J7 is the solid angle of the disk or ring as seen from the central star. The bolometric 

luminosity is the sum of the thermal luminosity and the scattered light luminosity of the dust 
grains. For dust grains with albedo to, the scattered light luminosity is ivL; the thermal luminosity 
is (1 — uj)L. These approximations assume a grey opacity and are reasonable for r < 1. For a disk 
or ring geometry, the solid angle is 0/47r = H/a, where H is the scale height and a is the radial 
distance from the central star. 

We calculate the radial optical depth through the grid to estimate the amount of absorption 
and scattering of starlight by dust particles in the grid. We divide the optical depth r through the 
grid into two pieces, measures the radial optical depth of bodies explicitly in the grid and 
measures the radial optical depth of smaller bodies. Wc adopt the geometric optics limit for bodies 
with X where A is the wavelength of observation. For the large bodies in the grid: 

N M 

Td = ^^Nikaik Aai (A24) 

i=l k=2 

where Nik is the number density of mass batch k in annulus i and CTik is the extinction cross-section. 
We adopt = ^Trr^^^ and a volume Vik = 47rajAaj-ffjfc. The optical depth is then independent of 
the width of the annulus: 

N M 2 

yy^^^kr^^ (A25) 

1=1 k=2 ' 

For the smaller bodies, we must make assumptions about the variation of n and H with 
particle size. In most calculations, the vertical scale height of the dust iJjfe is either constant or 
grows slowly for smaller particles. Adopting a constant Hik thus maximizes the optical depth of 
the smallest grains. The cumulative number density usually follows Nc oc with /? ~ 2.5 for 
calculations with fragmentation (Dohnanyi 1969; Williams & Wcthcrill 1994; Tanaka ct al. 1996). 
Gas and Poynting-Robertson drag preferentially remove the smallest particles from the grid; these 
processes should lower the exponent in the number density to /? « 1.5 when either pg or is large. 
Several test calculations confirm this dependence. We thus consider two cases, /3 = 2.5 - when 
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small particles dominate the opacity and /? = 1.5 - when large particles dominate the opacity. The 
optical depth is then 

N — oo 2 -'V 2 



j=i fe=i 

where 



1=1 k=l 1=1 



K = ^ = (A27) 

^ \ 500 {rmin/l l^m)-^/^ /? = 2.5 ^ ^ 

Radiation pressure limits the size of the smallest particle to radii ~ 1-10 /xm. We plan to address the 
large range in possible optical depths in future studies where we explicitly calculate the evolution 
of small particles. 

With these definitions for the optical depth, the bolometric luminosity from planetesimals in 
the grid is 

N M 2 

^^ = EE^- (A28) 

i=l k=2 ^ 

This expression is independent of the scale height H. The bolometric luminosity from small particles 
in the grid is 

in the low optical depth limit and 

TV 2 

Lmax = 500 (rmin/l ^imy'/^ ^ (A30) 



1=1 * 



in the high optical depth limit. 



Krivov et al. (2000) model the dynamical evolution of two small grain populations in debris 
disks. They derive size distributions for grains with rik < 1 mm for various assumptions for the 
collisional evolution of larger grains. Their calculations do not include velocity evolution. Krivov et 
al. (2000) show that collisions dominate Poynting-Robertson drag in debris disks with substantial 
optical depth, r > 10^^. If supported by calculations with velocity evolution, these results suggest 
that our maximum optical depth is appropriate for the early evolution of debris disk systems. 



A. 7. Tests of the Evolution Code 

Kenyon & Bromley (2001) describe various successful attempts to match results derived from 
the multi-annulus code with published calculations of velocity evolution. We have fully tested the 
new code in a variety of situations to verify that our algorithms conserve mass, angular momentum, 
and kinetic energy throughout an evolutionary sequence. The code reproduces previous results on 
particle growth in the Kuiper Belt (Kenyon & Luu 1999). Comparisons with previously published 
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calculations of terrestrial planet formation (e.g., Weidenschilling et al. 1997) will be described in 
a forthcoming publication. 
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Fig. 1. — Evolution of particle number (left panels) and eccentricity (right panels) for planetesimals 
in orbit around a 3 Mq star at 35 AU. The sixteen annuli in the grid for this ring model contain 
planetesimals with ro = 0.1-10 m, cq = 2 x 10~^, N{m) oc m~^, and a total mass of 0.667 
at t = (top two panels). The middle panels show the particle number and eccentricity at t 
= 1 Myr; the bottom panels show the particle number and eccentricity at t = 3 Myr. Despite 
considerable mass loss from shattering, coUisional damping reduces the orbital eccentricities to ~ 
10~^ on timescales of 1-10 Myr. Low velocity collisions then promote growth of larger objects 
through mergers. Style adapted from Weidenschilling et al. (1997). 
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Fig. 2. — Evolution of optical depth (upper panel) and reprocessed luminosity (lower panel) for the 
ring model of Figure 1. 




Fig. 3. — As in Figure 1, for planetesimals with Vmax = 10 km at i = (top panels), t = 1 Myr 
(middle panels), and t = 4 Myr (bottom panels. 
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Fig. 4. — Damping time (left panels) and dust luminosity (right panels) as a function of initial 
surface density for ring models at 35, 70, and 140 AU. Models with S = Sq have surface densities 
of solid material similar to a minimum mass solar nebula. All models have Vmin = 10 cm; model 
results for different values of Vmax are shown as indicated by the legend in each panel. 
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Fig. 5. — The ratio of final mass, Mf, to the initial mass, Mj, as a function of the initial eccentricity 
for ring models at 35 AU (open circles), 70 AU (filled circles), and 140 AU (open diamonds). Models 
with large cq lose most of their mass. 




Fig. 6. — As in Figure 1, for a ring model with gas drag. The initial gas to dust mass ratio is 10:1. 




Fig. 7. — As in Figure 6, for a ring model with an initial gas to dust mass ratio of 100:1. 
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Fig. 8. — Evolution of particle number (left panels) and eccentricity (right panels) for planetesimals 
in orbit around a 3 Mq star at 30-150 AU. The 64 annuli in the model grid contain planetesimals 
with ro = 0.1-10 m, cq = 2 x 10~^, N{m) oc to~^, and a total mass of 10 at i = (top two 
panels). The model grid extends from 30 AU at the inner edge to ~ 150 AU at the outer edge. 
Despite considerable mass loss from shattering, collisional damping reduces the orbital eccentricities 
to ~ 10~^ on timescales of 1-10 Myr. Low velocity collisions then promote growth of larger objects 
through mergers. 
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Fig. 9. — Evolution of the smallest and largest bodies in the full disk model of Figure 8. Top panel: 
eccentricity of the smallest bodies, = 10 cm, as a function of disk radius at t = 0, 1, 10, and 50 
Myr. The abrupt rise in e at the outer edge of the grid is a numerical artifact. Middle panel: scale 
height of bodies with = 10 cm as a function of disk radius at t = 0, 1, 10, and 50 Myr. Bottom 
panel: mass of the largest body in each annulus at i = 0, 1, 10, and 50 Myr. Large bodies at the 
inner edge of the grid are starting the runaway growth phase just as growth begins at the outer 
edge of the grid. 
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Fig. 10.— 



Evolution of dust luminosity for the full disk model of Figure 8. 
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Fig. 11. — Evolution of surface brightness for the full disk model of Figure 8. At t = 0, the surface 
brightness profile is a power law, / oc a~^, with p = 3.5. For t = 0-1 Myr, the surface brightness at 
the inner edge of the disk fades by an order of magnitude; the surface brightness at the outer edge, 
log ''i ^ 1.8-2.0, barely changes. The outer edge of the disk begins to fade at t = 1-10 Myr. For t 
= 10-50 Myr, the entire disk fades by an order of magnitude, and begins to produce planets. For 
t = 30 Myr and 50 Myr, the rise in surface brightness at the outer edge of the disk is a numerical 
artifact. 
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Fig. 12. — Model images for the full disk calculation of Figure 8. The panels show surface brightness 
distributions in the disk at selected times in the calculation. 



